BEC-BCS crossover and universal relations in unitary Fermi gases 
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The contact parameter in unitary Fermi Gases governs the short-range correlations, and high- 
momentum properties of the system. We perform accurate quantum Monte Carlo calculations 
with highly optimized trial functions to precisely determine this parameter at T=0, demonstrate 
its universal application to a variety of observables, and determine the regions of momentum and 
energy over which the leading short-range behavior is dominant. We derive Tan's expressions for 
the contact parameter using just the short-range behavior of the ground-state many-body wave 
function, and use this behavior to calculate the two-body distribution function, one-body density 
matrix, and the momentum distribution of unitary Fermi gases; providing a precise value of the 
contact parameter that can be compared to experiments. 
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The experimental realization and concurrent theoret- 
ical calculations of two- component unitary Fermi gases 
with short-range interactions offer a unique opportunity 
to test our understanding of strongly-interacting Fermi 
systems, and to study their structure and dynamics. The 
low-energy properties of the system are governed by the 
Bertsch parameter ^, the pairing gap i5, and have been 
studied extensively in the literature (TJ [2]. The short- 
range correlations of the system are, in contrast, governed 
by the many-body wave function at small interparticle 
separations, as encoded in the contact parameter C. 

Tan, in a series of papers [3], showed that in Fermi 
gases, if the effective range of the interaction is much 
smaller than any other length scale of the system, several 
universal relations occur, and related them to a param- 
eter C he called the contact parameter (see also Ref. [4 
for a review). In particular, the large momentum tail 
of the momentum distribution N{k), behaves like C/k'^. 
This same parameter gives the small distance behavior 
of the two-body distribution function, and the derivative 
of the ground-state energy with respect to the inverse of 
the scattering length. Having multiple phenomena that 
depend on a single universal parameter means that the 
parameter can be calculated and measured in multiple 
ways, and in the range of validity of the experiments 
and the calculations, they must give the same results for 
the parameters. We employ highly-optimized trial wave 
functions and accurate quantum Monte Carlo to calculate 
several of these observables, and to extract the contact 
parameter. We thus produce more accurate results for 
the leading behavior and simultaneously determine the 
regimes where it is dominant. 

Recently, several experimental measurements have 
measured C using a variety of techniques Values 
for the contact parameter have been calculated from pre- 
vious results from quantum Monte Carlo E] and other 
methods [7]. However, previous quantum Monte Carlo 
calculations give values of the contact at unitarity that 
disagree with each other at the 5 to 10 percent level. 



Here we show that if the calculations are carefully op- 
timized and extrapolated to zero range, our quantum 
Monte Carlo results agree with each other within statis- 
tical errors, less than 0.5 percent, giving clear numerical 
evidence of Tan's predicted universal contact parameter 
and its behavior around unitarity. Our results provide a 
benchmark prediction for low temperature experiments. 

We perform variational and fixed-node diffusion Monte 
Carlo (VMC and DMC) calculations of a system of a ho- 
mogeneous system of fermions interacting with a short 
range potential. Fixed-node diffusion Monte Carlo re- 
sults produce upper bounds to the ground-state energy 
of the system depending only upon the nodes (zeroes) of 
the trial wave function. We carefully optimize the trial 
wave functions, and obtain the best upper bounds to date 
for the ground state energy. Observables other than the 
ground-state energy are calculated by extrapolating the 
variational and mixed estimates: Ov = (5't|0|5't) and 
0,n = (*o|0|*t) , (O) « 20„i-0y. After suitable opti- 
mizations the extrapolations are very small. We calculate 
the contact parameter in several different ways and show 
they all give results consistent with each other and with 
recent experiments. 

Tan and later others [31 derived expressions for his 
contact parameter using a variety of methods. These 
results can be understood as coming from the behavior 
of the many-body wave function when two unlike spin 
particles are separated by a distance r small compared 
to the average particle separation, tq, but outside the 
range of the potential, R, 



(a is the two-body scattering length) which is Eq. 1 
in [3j and will be seen to be proportional to Tan's 
contact parameter C. The unlike spin two-body distri- 
bution function will be given by (r) in this same range 

guir) = A'{r-'-2a-'r-' + ...), (2) 
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where g^ir — oo) = ^ for an unpolarized system. The where Epc 



momentum distribution summed over both spins wiU also 
be dominated by this short range part of the wave func- 
tion, so for k much greater than the Fermi momentum 
kp but much less than the inverse potential range, we 
have 
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where dR indicates integration over ri, . . . , r^, and n is 
the number density. Fourier transforming the momen- 
tum distribution and the two-body distribution functions 
gives the behavior for the one-body density matrix (nor- 
malized to 1 at the origin) and the opposite spin static 
structure factor (which goes to | for k — > oo) of 

p(^)(r) = 1 -2TrnA'^r + . . . 
Snik) - o = r 
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Tan also related the contact parameter to the deriva- 
tive of the energy with respect to the inverse scattering 
length. Changing the scattering length by changing the 
potential with the mass fixed and using the Hellman- 
Feynman theorem [5] 



dE 
da-^ 



,3 / ^dv{r) 



(4) 



where E is the energy per particle. Since v(r) is 
nonzero only inside i?, where the two-body potential 
is very strong, g{r) can be replaced with /^(r) where 
^V^/(r) — v{r)f{r) and the integration taken over a 
sphere of radius R. Therefore 
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This only depends on /(r) around R, and using Eq. [T] 
the result is 
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The equation of state and therefore Tan's C [TUj are 
conventionally parametrized around unitarity as [5] 

E ^ C bv 
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j^p^ is the infinite system free gas en- 
ergy per particle. At unitarity we have several quantities 
related to C: 



p^^^r) ^ 1 - —Ckpr , N{k) 
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We use Quantum Monte Carlo (QMC) techniques to 
accurately solve the many-body ground state, and com- 
pute properties of the unitary Fermi gas. Our QMC cal- 
culations use the many-body Hamiltonian, 
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where the two-body interaction is a a short-range poten- 
tial taken only between opposite spin particles. At uni- 
tarity, vq — 1 and the effective range is r^- The scattering 
length and effective range can be tuned by changing vq 
and Tg. The limit of zero effective range (dilute system) 
is reached by taking r^, <^ rg, with vq = (3/(47rn))^/'^. 
The unitary limit is approached when vq <^ a where a 
is the scattering length of the two-body interaction. At 
unitarity the details of the interaction are not important, 
and the only scale of the system is given by its Fermi mo- 
mentum kp. The ansatz for the many-body trial wave 
function is the same as previously used in other QMC 
calculations [TT] : 

"^T = '[[fj{nj') $BCS, $BCS =-4[0(riiO0(r22')--0(?-nn 

ij 

where A antisymmetrizes the like spins, and the un- 
primed coordinates are for up spins and the primed are 
for down spins and n = N/2. The pairing function is 

(/)(r) ^ I3{r) + a{kl) exp[ik„ • r] , 

n 

^{r)^P{r)+PiL~r)-2f3{L/2) , 

I3{r) = [1 + cbr] [1 - exp(-d6r)] "'""'"['^'"^ . (9) 

dor 

The function f3{r) has a range of L/2, the value of c is 
chosen such that it has zero slope at the origin. 

The variational wave function has been carefully op- 
timized; in particular we optimize the pairing orbitals 
entering in the wave function by using VMC to minimize 
the energy p2]. The fixed-node DMC energies do not 
depend on the Jastrow function fj. Our simulations are 
performed with 66 particles in a periodic box, and we 
study the effect of the effective range of the interaction 
by changing r^ and extrapolating to the r^ — > limit. 
The results of 66 particles is very close to the infinite 
limit [13]. Careful optimization of the variational wave 
function significantly improves the energy upper bounds. 
At unitarity, the best previous QMC results using 66 
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FIG. 1. (color online) Energy per particle in the BEC-BCS 
crossover regime in units of Efg as a function of the scattering 
length a. The QMC points are the results of extrapolations 
to re — > limit. In the inset we show the extrapolation at 
unitarity. 
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FIG. 2. (color online) The calculated spin-opposite two-body 
distribution function at small distance r; the effective range 
of the interaction is kp ~ 0.08 (vertical black dashed line). 
The VMC (red), mixed (green) and extrapolated (blue) re- 
sults are shown. The extrapolated QMC results are used to 
fit the function giving b = 1.2678. In the inset we show the 
same functions multiplied by (kpr)^ ■ 



particles are ^ = 0.42(1) fixing rg/ro ~ 0.08 ^j, and 
f = 0.42(1) using r^/rQ « 0.01 [14 . Our new estimate 
is C = 0.4069(5) and ^ = 0.3923(4) with r^/ro ~ 0.07 
and 0.02 respectively. The parameters for at unitarity 
are b = 0.5kp, d = 5 and the nonzero a{k'^) are given in 
Table |l] Improved optimization of the trial wave function 



lowers the fixed-node energy by 4-7%. Careful extrap- 
olation to Te — limit is also important. We show an 
example at unitarity in the inset of Fig. [l] where we plot 
QMC points at different effective ranges, and their ex- 
trapolation. Using more points and a more complex fit 
typically provides a somewhat lower upper bound to the 
energy such a correction is about 0.002 to ^. 

We optimized the many-body wave function for sys- 
tems with different scattering lengths and for each value 
of kpa we repeated the extrapolation of rg. Our re- 
sults of S,{kpa) are shown in Fig. [T] Fitting the QMC 
points shown in Fig. [l] gives the values f — 0.383(1), 
C = 0.901(2) and ly = 0.49(2). Using Eq. [5]we predict 
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An alternative direct method for calculating the con- 
tact can be obtained by computing correlation functions 
at unitarity. For example, the pair distribution function 
is shown in Fig. [2] where we compare the VMC result 
with the mixed estimate computed with DMC. The two 
results are almost identical and differences appear only 
for very small distances. The value of C is obtained by 
fitting g^i{kpr) in the range <C r ^ kp^ using the 
function a + b/r"^. The fit gives b = 1.2678(1). Using 
Eqs. [2] and [sj gives the value for C = 0.897(2) in good 
agreement with the result extracted from Eq. [6] 

The calculated radial one-body density matrix 
p^^^kpT) is shown in Fig. [3] using VMC and the mixed 
DMC results. Again the results are nearly identical, with 
strikingly linear behavior over a large range of small kpr 
values. The fit gives C = 0.895(16) again in good agree- 
ment with the equation of state result. The calculated 
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0.00198 

1 0.00250 

2 0.00194 

3 0.00081 

4 0.00033 


5 0.000190 

6 0.000200 

8 0.000167 

9 0.000163 

10 0.000120 



TABLE I. The optimized plane wave coefficients at unitarity 
for the pairing function. 



FIG. 3. (color online) The radial one-body density matrix, 
symbols and kp « 0.08 as in Fig. [2] A line showing the 
linear fit with c — 0.2685 is also shown, the dominant short- 
range behavior is accurate up to approximately kpr ~ 3. 

momentum distribution is shown in Fig. |4] The momen- 
tum distribution and the one-body density matrix are 
each other's Fourier transform. The only difference in 
our calculations are that the angular average has been 
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should be useful to future experiment in extracting the 
contact behavior and leading corrections. 
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FIG. 4. (color online) The calculated momentum distribution 
summed over both spins multiplied b y fc^ /fep showing the fc"* 
tail. Dashed line show 2C/kp of Eq. 
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done in real space for the one-body density matrix to 
give the radial one-body density matrix, while the mo- 
mentum distribution is calculated for the k vectors that 
correspond to the periodic simulation cell. The extrac- 
tion of the fc^ tail is rather noisy; using the radial one- 
body density matrix gives a more accurate fit. From our 
results it appears that the contact term dominates the 
behavior for k > 2k p. Our asymptote is consistent with 
the value 0.229(1) expected from C = 0.901(2) (dashed 
line in Fig. [4]). 

Recent experiments have measured the contact pa- 
rameter from the equation of state .5;, momentum dis- 
tribution directly using ballistic expansion and indi- 
rectly through the rf line shape and photoemission spec- 
troscopy [S], and from the static structure factor [S]. 
Navon et al. [5] extracted a value of C = 0.93(5) from 
their equation of state measurements. Our best value 
of C = 0.901(3) is well within their experimental errors. 
Kunhle et al. calculate a slope of S{k) versus kp/k 
at large k for l/ikpa) = of 0.75(3) at T = 0.10(2)rF, 
giving a value of C = 0.80(3), while Stewart et al. give 
values somewhat away from unitarity which also give C 
lower than our value. 

In conclusion, wc have used Quantum Monte Carlo 
techniques to study the short-range correlations of uni- 
tary Fermi gases as encoded in Tan's contact parame- 
ter. The extractions from various observables all give the 
same result within statistical errors. These Monte Carlo 
methods give particularly low variance values for the en- 
ergy of the system and with minimal bias. Therefore 
extracting the contact parameter from the equation of 
state is the simplest and most reliable. However, we have 
shown that its value extracted from the two-body radial 
distribution function, the one-body radial density matrix, 
and the momentum distribution also give the same re- 
sults albeit with somewhat larger error bars. For each of 
these quantities we have also determined the regime over 
which the leading contact behavior is dominant, which 



[1] S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. 

_Phys. 80, 1215 (2008); J. Carlson and S. Reddy, Phys!" 

Rev. Lett. 95, 060401 (2005); Phys. Rev. Lett. 100, 
150403 (2008). 

[2] A. Schirotzek, Y.-i. Shin, C. H. Schunck, and W. Ketterle, 
Phys. Rev. Lett. 101, 140403 (2008); D. Lee, Phys. Rev. 
C 78, 024001 (2008). 

[3] S. Tan, Ann. Phys. 323, 2952 (2008); 323, 2971 (2008); 
323, 2987 (2008). 

[4] E. Braaten(2010), arXiv: 1008.2922 

[5] N. Navon, S. Nascimbene, F. Chevy, and C. Salomon, 

Science 328, 729 (2010); E. D. Kuhnle, H. Hu, X.-J. 

Liu, P. Dyke, M. Mar k, P. D. Drum mond, P. Hannaford, 

and C. J. Vale, Phys. Rev . Lett7| 105, 070402 (2010); 

J. T. Stewart, J. P. Gaeblef7 TT E. Drake, and D. S. 

Jin, 104, 235301 (2010); H. Hu, X.-J. Liu, and P. D. 

Drummond(2010), ar Xiv:1011.3 845 
[6] C. Lobo, I. Carusotto, S. Giorgini, A. Recati, and 

S. Stringari, Phys. Rev Lett.| 97, 100405 (2006); 



R. Combescot, S. Giorgini, and S. Stringari, Europhys. 
Lett. 75, 695 (2006); J. E. Drut, T. A. Lahde, and T. Ten 
arXiv: 1012.5474; S.-Q. Su, D. E. Sheehy, J. Moreno, and 
M. Jarrell, Phys. Rev. A 81, 051604 (May 2010); T. Abe 
and R. Seki, Phys. Rev. C 79, 054003 (2009). 
[7] F. Palestini, A. Perah, P. Fieri, and G. C. Strinati, Phys. 
Rev. A 82, 021605 (2010); J.-W. Chen and E. Nakano, 
75, 043620 (2007); D. Lee, Eur. Phys. J. A 35, 171 
(2008). 

[8] E. Braaten and L. Platter, Phys. Rev. Lett. 100, 205301 
(2008); R. Combescot, F. Alzetto, and X. Leyronas, 
Phys. Rev. A 79, 053640 (2009); F. Werner, L. TarrueU, 
and Y. Castin, Eur. Phys. J. B 68, 401 (2009). 

[9] H. Hellmann, Einfiihrung in die Quantenchemie 
(Leipzig: Franz Deuticke, 1937) p. 285; R. P. Feynman, 
Phys. Rev. 56, 340 (1939). 
[10] Some authors define the contact as an extensive quantity 
C = fiC, where f2 is the volume, and report the unitless 
intensive quantity = 3n^ p- . 

[11] J. Carlson, S.-Y. Chang, V. R. Pandharipande, and K. E. 
Schmidt, Phys. Rev. Lett. 91, 050401 (2003); S. Y. 
Chang, V. R. Pandh aripande, J. Carlson, and K. E. 
Schmidt, 'Phys. Rev. A 70, 043602 (2004); A. Gezerlis, 
S. Gandolfi, K. E. Schmidt, and J. Carlson, Phys. RevT] 



Lett. 103, 060403 (2009). 
[12] S. Sorella, Phys. Rev. B 64, 024512 (2001). 
[13] M. McNeil Forbes, S. Gandolfi, and A. Gezerhs(2010), 

,arXiv:1011.2197. 



5 



[14] G. E. Astrakharchik, J. Boronat, J. CasuUeras, and 
S. Giorgini, Phys. Rev. Lett. 93, 200404 (2004). 



